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The dynamics of a phase-separated two-component Bose-Einstein condensate are investigated, in 
which a bubble of one component moves through the other component. Numerical simulations of 
the Gross-Pitaevskii equation reveal a variety of dynamics associated with the creation of quantized 
vortices. In two dimensions, a circular bubble deforms into an ellipse and splits into fragments with 
vortices, which undergo the Magnus effect. The Benard-von Karman vortex street is also generated. 
In three dimensions, a spherical bubble deforms into toruses with vortex rings. When two rings are 
formed, they exhibit leapfrogging dynamics. 

PACS numbers: 67.85.Fg, 67.85.De, 47.32.cf, 47.32.ck 



I. INTRODUCTION 

A droplet of ink falling into water deforms from a 
sphere to a toroidal shape due to the formation of a 
vortex ring p], Q . An air bubble rising in water ex- 
hibits complicated dynamics such as zigzag and spiral 
motion 0, Q . Such phenomena are caused by the inter- 
action between the droplet (ink) or the bubble (air) with 
the medium (water), where the former moves through the 
latter. The subject of the present paper is the dynamics 
of phase-separated two-component superfluids in similar 
situations, that is, a bubble of one component moving 
through the other component. 

A variety of dynamical properties have been observed 
in two-component Bose-Einstein condensates (BEC) 
Two immiscible components prepared in a spatially over- 
lapped state separate dynamically [6| and develop into a 
domain structure [7]. A similar structure is observed in 
a miscible system with counterflow [8]. In rotating two- 
component BECs, a square vortex lattice has been ob- 
served [9] . These dynamical properties of two-component 
BECs have been studied theoretically p~ol-fl5| . Theoreti- 
cal predictions have been made on various structures in 
two-component BECs; for example, Skyrmions [16|, soli- 
tary wave complexes [TtJ , and vortex bright solitons [l8| . 
Various interface instabilities known in classical fluid me- 
chanics are predicted to emerge in two-component BECs 
with an interface, namely the Rayleigh-Taylor instabil- 
ity [H, Hj, the Kelvin-Helmholtz instability [H El, 
and the Richtmyer-Meshkov instability [23j. The be- 
havior of two-component BECs strongly depends on the 
miscibility, which is determined by the intra- and inter- 
component interactions. Recently, these interactions in 
two-component systems have been controlled using the 
Feshbach resonance [24[, and phase separation dynamics 
have been studied in a controlled manner [25|, [26| . 

In the present paper, we investigate the dynamics of a 
bubble in a phase-separated two-component BEC, where 
a small fraction of one component moves through the 
other component. We show that the system exhibits a 
rich variety of dynamics depending on the parameters 
and dimensionality. In two dimensions (2D), a circular 
bubble at rest deforms into an ellipse as it accelerates and 



breaks into pieces with quantized vortices. For strong 
phase separation, a vortex street is formed in the wake of 
the moving bubble, as in the Benard-von Karman vortex 
street in a single-component BEC [27]. In three dimen- 
sions (3D), a spherical bubble deforms into a torus, as 
in classical fluids p], Q . The significant difference from 
classical fluids is that the toroidal bubble is accompanied 
by a quantized vortex ring. When two or more rings are 
generated, they leapfrog each other. 

This paper is organized as follows. Section HI1 provides 
a formulation of the problem. Section fill Al numerically 
demonstrates the dynamics of a 2D system and Sec. lIIIBl 
analyzes the deformation of a bubble. Section IIVI per- 
forms full 3D numerical simulations. Section [V] gives 
conclusions to this study. 



II. FORMULATION OF THE PROBLEM 

We study the dynamics of a two-component BEC using 
the zero-temperature mean-field theory. The dynamics of 
the system are described by the Gross-Pitaevskii (GP) 
equation given by 
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where ipj is the macroscopic wave function, rnj is the 
atomic mass, and Vj is the external potential for the jth 
component (j = 1,2). The interaction parameters g^> 
are defined as 



gjji = 27rh 2 a jj f(m j 1 + m-, 1 ), 
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where a^y is the s-wave scattering length between the 
atoms in the jth and j'th components. 

We assume that the interaction parameters satisfy 
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and therefore the two components are phase sepa- 
rated [28]. We also assume that the external potentials 
for both components are initially absent, V\ = V2 = 0. 
We consider an initial state in which a small fraction of 
component 2 is localized and surrounded by component 
1, and ipi is uniform far from component 2. The ground 
state is therefore a circular (2D) or spherical (3D) "bub- 
ble" of component 2 located in the sea of component 1. 
At t = 0, we start to exert a force F on the bubble by 
applying a potential V2 = Fx to component 2. The bub- 
ble then moves in the — x direction. This is done, e.g., by 
applying a magnetic field gradient to a system in which 
the atoms in component 2 have a magnetic dipole mo- 
ment in the direction of the magnetic field while those in 
component 1 do not. 

We numerically solve Eq. (pQ) by the pseudo-spectral 
method. The initial state is the ground state for V\ = 
V2 = obtained by the imaginary-time propagation 
method. We add a small amount of white noise to the ini- 
tial state to break the numerically exact symmetry. The 
boundary in the numerical calculation is set far from the 
relevant region so that the periodic boundary condition 
does not affect the results. In the following analysis, we 
assume mi = rri2 = m and gu = #22 = g to reduce the 
number of parameters. 

III. TWO-DIMENSIONAL SYSTEM 
A. Dynamics 

We first demonstrate the numerical results for the 2D 
system. The wave function is assumed to have the form 
ifjj(r,t) = ifjj±(x,y,t)(j)(z), where the dynamics of the 
normalized wave function cj)(z) are frozen. Integrating the 
GP equation with respect to z, the effective interaction 
strength is found to be g^y J \ (j)\ A dz. Hence, the healing 
length is ^2D = ft/(mgri2B J \(\>\^dz) x l 2 and the sound ve- 
locity is Vg D = (gn2B J \ (})\ A dz /m) 1 / 2 : , where U2B is the 2D 
density |^i±| 2 f ar from the bubbles. Normalizing length 
and time by and ^d/v s 2D , the intra-component in- 
teraction parameter in the GP equation becomes unity 
and the relevant interaction parameter is only gvijg- 

Figure [1] shows the time evolution of the density dis- 
tribution. The initial state is a circular bubble as shown 
in Fig. [2 (a). At t = 0, we apply a potential V2 = Fx 
to component 2, and the bubble is accelerated in the — x 
direction. We find that the bubble is deformed in an el- 
liptic shape as shown in Fig. [I] (b). After the potential 
V 2 is switched off at i = fr^ D /f 2D = 200, the bubble 
moves at a constant velocity keeping the elliptic shape 
as shown in Fig. [1] (c). Such elliptic deformation of a 
bubble is known in classical fluid mechanics [3] and will 
be analyzed in Sec. IIII Bl 

Figures [2] (a)-[2] (c) show the dynamics for the same 
parameters as those in Fig. [TJ where the force is exerted 
for t > 0. The bubble is first deformed elliptically as in 
Fig. [D (b) and then splits into two bubbles [Figs. [2] (b) 




FIG. 1: (Color online) Time evolution of the normalized 
density distribution hj = | 2 / n 2D for gi2/g = 1.25. From 
t = t^ D /f 2 D = to t = 200, the force F = F& B /(hv* D ) = 
1.58 x 10 -3 is exerted on component 2 in the — x direction. 
The amount of component 2 is J \ip2\ 2 dxdy /(tizd^zd) — 205. 
The field of view is 158 x 126 in units of ^2D. 

and [2](c)]. We note that component 1 contains singly- 
quantized vortices with opposite circulations at the split 
bubbles as indicated by the arrows in Fig. [2] (c). After 
the split, the two bubbles move away from each other as 
shown by the dashed lines in Fig. [2] (d), and eventually 
move in the directions perpendicular to the force. This 
behavior is due to the Magnus force on the vortices given 
by 

F M = pK x v, (4) 

where p is the mass density, n is the vorticity, and v 
is the velocity of the vortex. In this case, p ~ mri2D 
and k = 27rh/m. The velocity of the bubbles in the ±y 
directions is estimated by equating the Magnus force in 
the +x direction and the external force on the bubble in 
the —x direction as 

mn2Y^^ L v y =F I \ijj2±\ 2 dxdy, (5) 

TTl J 

where the integration is taken within one of the bubbles. 
For the parameters in Fig. [2j the velocity is estimated 
from Eq. (|5j) to be v y /v^ D ~ 0.026, which is in good 
agreement with the velocity ~ 0.026 along the dashed 
line in Fig. [2] (d). If the external potential for compo- 
nent 2 is switched off after the bubble splits, the two 
bubbles move in the — x direction at a constant velocity 
~ ft/(27r<i) with a constant distance d maintained be- 
tween the bubbles, as shown by the dotted lines in Fig. [2] 
(d). Such a two-component vortex is suitable for study- 
ing the Magnus effect, since we can exert a force on the 
vortex in a controlled manner. 

Figures [3] ( a) -[3] (f) shows the time evolution of the den- 
sity profile of component 1, where the force is stronger 
than that in Figs. □ and [2 At t = 300, the bubble re- 
leases two fragments, at which point component 1 has 
a vortex- ant i vortex pair [two white circles in Fig. [3] (a)]. 
The preceding bubble then splits into two, at which point 
a vortex- ant ivortex pair is also contained in component 
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FIG. 2: (Color online) (a)-(c) Time evolution of the density 
hj and phase <j>j = arg^j profiles. The force F = 1.58 x 
10 -3 is kept constant for t > 0. Other parameters are the 
same as those in Fig. [1] The arrows indicate the directions 
of circulation of the vortices, (d) Stroboscopic trajectories of 
component 2, where the images are taken at an interval of 
5t — 1000. The force is constant for the images lying on the 
dashed lines, and the force is switched off at t — 1000 for the 
images lying on the dotted lines. The field of view is 158 x 126 
in (a)-(c) and 316 x 253 in (d) in units of ^2D- 



(d) 7 =2600 



7=4500 









11.2 


<2) 
o 




© • 






© 




® 










1 



(g)7=500 

FIG. 3: (Color online) (a)-(f) Time evolution of the density 
profile of component 1, where the force F — 4.74 x 10 ~ 3 is 
exerted on component 2 from t = to t = 500. The other 
conditions are the same as those in Fig. [T] The two vortex 
pairs are marked by black and white circles to follow their 
motion, (g) Phase profile of component 1 at t = 500. The 
field of view in (a)-(g) is 316 x 126 in units of <^2D- 



bubble changes in time. Occasionally, vortices are shed 
in pairs periodically as shown in Fig.|4](b), which is rem- 
iniscent of the Benard-von Karman vortex street in a 
single-component BEC [27j . Unlike vortex shedding in a 
single-component BEC, the bubble gradually diminishes 
because the cores of the released vortices are occupied by 
component 2 as shown in the bottom panel of Fig. |4] (b). 
The condition thus continuously changes and the vortex 
street generation does not last long, since the parameter 
region for the vortex street generation is narrow [27| . 
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(a) t = 300 



(b) t = 500 




(e) 7= 3600 



1 [black circles in Fig.[3](b)]. After that, the pair marked 
by the white circles passes through that marked by the 
black circles, and the two pairs leapfrog each other back 
and forth as shown in Figs. [3] (c)-|3] (f). 

Figure [H shows the cases for strong force and strong 
phase separation. Since the interface tension is large for 
strong phase separation [29|, [30| , the bubble is hard to 
split as in Fig. [2j Consequently, the bubble is acceler- 
ated beyond the critical velocity for vortex creation and 
the vortices are shed in the wake as shown in Fig. HJ 
Figure [4] (a) shows a typical behavior of the bubble, in 
which the Magnus force bends the trajectory of the bub- 
ble while it contains vortices and then the bubble moves 
in a winding fashion since the vorticity contained in the 



B. Analysis of bubble deformation 

To understand the elliptic deformation of the 2D bub- 
ble shown in Fig. [TJ we perform a simple analysis. 

In this subsection, we assume that the thickness of 
the interface between the two components is negligi- 
ble due to the strong phase separation, and the effect 
of the interface is expressed only by an interface ten- 
sion coefficient a. The wave function of the stationary 
state of component 1 outside the bubble is written as 

1/2 

= n i ( r ) exp[i#i(r) — ifiit/h], where we con- 
sider the problem in the frame moving with the bubble. 
Substituting this expression into the GP equation (pQ), we 
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FIG. 4: (Color online) Snapshots of the density and phase 
profiles for gw/g = 3.3. The force on component 2 is F = 
6.96 x 10" 3 in (a) and F = 6.01 x 10" 3 in (b). The black 
crosses in (a) indicate the positions of the bubble at an interval 
of At = 200. The field of view is -869£ 2 d < x < -474£ 2 d and 
-158^2D < y < 158£ 2 d in (a) and -869£ 2 d < x < -474£ 2 d 
and — 63^2D < V < 63£ 2 d in (b), where the initial spherical 
bubble is located at x = y = 0. The amount of component 2 
is the same as that in Fig. [1] 



obtain 



h 2 v 2 v / M^) | i 
2m V^W 2 



V - Mr>i(r)] = 0, (6) 
mv 2 (r) + tfiini (r) = m, (7) 



where v\ = TiV0i/m. We also assume that the system 



is almost incompressible, i.e., rij(r) 



Srij(r) with 



Srij(r) <C rij. Equations (j6j) and ([7]) then become 

V-«i(r) - 0, (8) 



Pi( r ) H — nimv 2 (r) ~ const., 



(9) 



where Pi(r) — gun 2 (r)/2 is the pressure. Equation ([9]) 
corresponds to Bernoulli's equation in classical fluid me- 
chanics. 

The shape of the bubble is approximated to be an el- 
lipse given by x 2 /a 2 + y 2 /b 2 = 1, and we estimate a and 
b considering the pressure inside and outside the bub- 
ble at r a = (a, 0) and = (6,0). Since v\ = at the 
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FIG. 5: (Color online) Oblateness ft of an elliptic bubble as 
a function of the velocity vq. The solid curve shows Eq. ([T2]) 
with the interface tension coefficient a in Eq. JT3]). The dotted 
curve shows Eq. (|12|) with the left-hand side multiplied by 1.2. 
The circles are numerically obtained from the GP equation. 
The inset shows an example of the density distribution of 
component 2 for a stationary state. 



stagnation point r a , Eq. (j9]) gives 

1 



Pi(r a ) = Pi(r b ) 



-nimv^rb) 



(10) 



For inviscid and incompressible flow, the velocity on ei- 
ther side of an elliptic obstacle is v\{rb) = vo(a+b)/a (3lj . 
where is the velocity at infinity. Applying Laplace's 
formula [32[ to the inner pressure at r a and r&, we have 

Pi(r a ) + <*^ = Pi(n) + a\, (11) 
b z a z 

where a/b 2 and b/a 2 are the curvatures of an ellipse at r a 
and r&, respectively. Using Eqs. ([TU]) and ([IT]) , we obtain 
an expression that estimates the oblateness of the bubble 

as 



V/ ^(l + /3) 2 



(12) 



where A = nab is the area of the ellipse and j3 = a/b. 
The left-hand side of Eq. ([T2]) corresponds to the Weber 
number in fluid mechanics. 

We compare Eq. ([T2]) with the numerical calculation. 
We numerically solve Eq. (pp) with V\ = V2 = and 
—ivod x ipj added to the right-hand side. The imaginary 
time propagation of this equation gives the stationary 
state in the frame moving in the — x direction with ve- 
locity vq. The circles in Fig. [5] plot the ratio of the size of 
the bubble in the minor axis to that in the major axis as 
a function of the velocity. The solid curve in Fig. [5] shows 
Eq. ([T2]) with the interface tension coefficient [29|, [30| , 
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FIG. 6: (Color online) (a)-(d) Time evolution of the iso- 
density surface of component 2 in the 3D simulation. The 
size of the frame is 178 x 63 x 63 in units of £. (e) Den- 
sity hi = |^i| 2 /^3D and phase <\>\ — arg^i profiles of com- 
ponent 1 on the z — cross section at t = tv s /^ = 560. 
The field of view is 178 x 63 in units of £. The param- 
eters are g 12 /g = 1.1, F = F£ 2 /(hv B ) = 3.2 x 10" 3 , and 
/|V>2| 2 dr/(n 3 D£ 3 ) = 3.2 x 10 3 . 



We find that the solid curve deviates from the circles. 
The deviation of the analytic result from the numerical 
result may be due to the approximations and assump- 
tions made in the analysis. This deviation can be com- 
pensated if we modify Eq. ([T2]) . The dotted curve in 
Fig. [5] shows Eq. (fl2]) with the left-hand side multiplied 
by 1.2, which is in good agreement with the circles in 
Fig. El 



IV. THREE DIMENSIONAL SYSTEM 

In the numerical simulations in 3D, we normalize 
length and time by £ = %/ '{mgn^u) 1 ' 2 and £/v S: where 
usd is the atom density far from the bubbles and v s = 
(gnsn/m) 1 / 2 is the sound velocity. 

Figures [6] (a)-|6] (d) show typical dynamics of a bub- 
ble in 3D. The initial state is a spherical bubble [Fig. [6] 
(a)], which deforms into a ellipsoidal shape as it is ac- 
celerated [Fig. [6] (b)]. The bubble then takes a toroidal 
shape [Fig. [6] (c)], whose radius increases in time [Fig. [6] 
(d)]. Since the amount of component 2, f\ip2\ 2 dr, is 
conserved, the torus becomes thin as its radius increases. 
Figure [6] (e) shows the density and phase profiles on the 




FIG. 7: (Color online) Time evolution of the isodensity 
surface of component 2 in the 3D simulation. The force is 
F = 9.5 x 10 ~ 3 and the other parameters are the same as 
those in Fig. [6] The two objects are color coded for clarity. 
The size of the frame is 221 x 63 x 63 in units of £. 



cross section of z = at t — tv s /£ = 560. We find that 
component 1 contains a quantized vortex ring along the 
torus of component 2. If the force on component 2 is 
switched off after the ring is formed, the ring propagates 
in the — x direction at a constant velocity without ex- 
pansion of the radius. The velocity roughly agrees with 
v ~Ti/ (2mR) \n(SR/a) [33[, where R is the radius of the 
torus and a is the radius of the vortex core. 

In contrast to a vortex ring in a single component su- 
perfluid [34|, we can manipulate the vortex ring by ap- 
plying a potential to the torus of component 2 using a 
magnetic field gradient or laser field. We also note that 
a Skyrmion [16[ can be created if we imprint the phase 
arg(z ± iy) on the torus of component 2 using, e.g.^the 
Raman transition with Laguerre-Gaussian beams [35| . 

When the force F is strong, we can create multiple 
rings as shown in Fig. [7J where the two rings are color 
coded to distinguish them. The initial state is the same 
as that in Fig. [6] (a). The bubble first releases a ring 
backwards [Fig. (a)] and then the front bubble also 
deforms into a ring [Fig. [71 (b)], resulting in double vor- 
tex rings with the same circulation. After that, the rear 
ring passes through the front ring [Figs. 13(c) and [3(d)], 
and this overtaking is repeated [Figs. [3 (d) and [3 (f)], 
which is the 3D version of the behavior in Fig. [3l Such a 
leapfrogging behavior of vortex rings was first predicted 
inRef.0. 



V. CONCLUSIONS 

We have investigated the dynamics of phase-separated 
two-component BECs, in which a "bubble" of one com- 
ponent propagates in the other component. We studied 
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2D systems in Sec. IIIII When the velocity of the bub- 
ble is small, the circular bubble deforms into an elliptic 
shape, which travels with a constant velocity if the force 
is switched off (Fig. [1]). The elliptic deformation of the 
bubble was analyzed in Sec. IIII B[ For a large velocity, 
the bubble of component 2 splits into two or more frag- 
ments, where component 1 contains quantized vortices 
(Figs. [2] and [3]). The trajectories of bubbles containing 
vortices are then affected by the Magnus force [Fig. [2] 
(d)]. When the force and phase separation are strong, 
vortex shedding occurs instead of the split. The bubble 
drifts due to the Magnus force and sometimes generates 
the Benard-von Karman vortex street (Fig. |4]). For the 
3D system studied in Sec. [Vj we found that the spheri- 
cal bubble changes to a toroidal shape, where component 
1 contains a quantized vortex ring (Fig. [6j). When two 
vortex pairs or two vortex rings are created, they exhibit 
a leapfrog behavior (Figs. [3] and [7)). 



We have thus shown that bubbles in two-component 
BECs exhibit a rich variety of phenomena. When the 
bubble of component 2 splits in 2D or becomes a torus 
in 3D, quantized vortices are generated in component 
1, whose cores are occupied by component 2. We can 
therefore exert a force on the vortices in component 1 in 
a controlled manner through the force on component 2, 
which may allow manipulation of vortices in a BEC. 
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